{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "baba5557",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.stats import norm\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e55631a",
   "metadata": {},
   "source": [
    "# Identification of the parameters $\\{p, \\phi, \\alpha_0, \\alpha_1\\}$\n",
    "\n",
    "The model is\n",
    "$$\n",
    "\\max_{m_i} \\frac{(y_i-pm_i)^{1-\\sigma}}{1-\\sigma} + \\phi (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "The FOC is \n",
    "$$\n",
    "p (y_i-pm_i)^{-\\sigma} = \\phi\\alpha_0 \\exp(-\\alpha_1 - \\alpha_0 m_i)\n",
    "$$\n",
    "This gives a solution $m_i = {\\cal M}(y_i,p|\\sigma, \\phi, \\alpha_0, \\alpha1)$\n",
    "\n",
    "\n",
    "### US\n",
    "For the US, we have $p=1$ and $\\sum_i \\omega_i y_i \\equiv y = 1$. Therefore, the moments restrictions  \n",
    "$$\n",
    "\\Psi_1 = \\frac{p \\sum_i \\omega_i m_i}{\\sum_i \\omega_i y_i}\n",
    "$$\n",
    "$$\n",
    "\\Psi_2 = \\sum_i (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "$$\n",
    "\\Psi_3 = \\frac{(1-\\exp(-\\alpha_1 - \\alpha_0 m_{\\max}))}{(1-\\exp(-\\alpha_1 - \\alpha_0 m_{\\min}))}\n",
    "$$\n",
    "allows to identify $\\{\\phi,\\alpha_0,\\alpha_1\\}$\n",
    "\n",
    "\n",
    "### Europe\n",
    "\n",
    "For Europe, we keep $\\{\\phi,\\alpha_0,\\alpha_1\\}$ found in the previous step, and the moments restrictions\n",
    "$$\n",
    "\\Psi_1 = \\frac{p \\sum_i \\omega_i m_i}{\\sum_i \\omega_i y_i}\n",
    "$$\n",
    "$$\n",
    "\\Psi_2 = \\sum_i (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "allows to identify $\\{p,\\alpha_1\\}$, given that a calibration of $\\sum_i \\omega_i y_i \\equiv y$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e3d32d8",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5f76a0af",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Function for finding the decision rule for m\n",
    "\n",
    "def medexp(m, *param_m):\n",
    "    p, y, sigma, alpha0, phi, alpha1 = param_m \n",
    "    return p*( (y - p*max(0,min(m,y/p)))**(-sigma) ) - alpha0*phi*np.exp(-( alpha1 + alpha0*max(0,min(m,y/p)) )) \n",
    "\n",
    "# Function for finding parameters (alpha0, phi, alpha1) using US moments' restrictions\n",
    "def hetero_calib(pd_c, *param_c):\n",
    "    \n",
    "    Psi1, Psi2, Psi3, p, yy, sigma = param_c \n",
    "    \n",
    "    msolv = np.zeros(yy.size)\n",
    "    \n",
    "    for ii in range(yy.size):\n",
    "        param_cm  = (p, yy[ii], sigma, pd_c[0], pd_c[1], pd_c[2])\n",
    "        m0        = .1\n",
    "        msolv[ii] = max(0,fsolve(medexp, m0, args=param_cm))\n",
    "    \n",
    "    return [Psi1 - np.mean(msolv), \n",
    "            Psi2 - np.mean(1-np.exp(-pd_c[2]-pd_c[0]*msolv)),\n",
    "            Psi3 - (1-np.exp(-pd_c[2]-pd_c[0]*msolv[-1]))/(1-np.exp(-pd_c[2]-pd_c[0]*msolv[0])) ] \n",
    "\n",
    "# Function for finding (p, alpha1) using European moments' restrictions\n",
    "def hetero_calib_eu(pd_s, *param_s):\n",
    "    \n",
    "    Psi1, Psi2, yy, sigma, phi, alpha0 = param_s \n",
    "    \n",
    "    msolv = np.zeros(yy.size)\n",
    "    \n",
    "    for ii in range(yy.size):\n",
    "        param_sm  = (pd_s[0], yy[ii], sigma, alpha0, phi, pd_s[1])\n",
    "        m0        = .1\n",
    "        msolv[ii] = max(0,fsolve(medexp, m0, args=param_sm))\n",
    "    \n",
    "    return [Psi1 - (pd_s[0]*np.mean(msolv))/(np.mean(yy)), \n",
    "            Psi2 - np.mean(1-np.exp(-pd_s[1]-alpha0*msolv))] \n",
    "\n",
    "# Function for finding the price, given all other parameters (alpha0, phi, alpha1)\n",
    "def hetero_new_price(p_new, *param_s):\n",
    "    \n",
    "    Psi1, yy, sigma, phi, alpha1, alpha0 = param_s \n",
    "    \n",
    "    msolv = np.zeros(yy.size)\n",
    "    \n",
    "    for ii in range(yy.size):\n",
    "        param_sn  = (p_new, yy[ii], sigma, alpha0, phi, alpha1)\n",
    "        m0        = .05\n",
    "        msolv[ii] = max(0,fsolve(medexp, m0, args=param_sn))\n",
    "    \n",
    "    return Psi1 - (p_new*np.mean(msolv))/(np.mean(yy)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8bf25d4",
   "metadata": {},
   "source": [
    "## Calibrated parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "49f43316",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma  = 2.0\n",
    "Ncount = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf29370a",
   "metadata": {},
   "source": [
    "## Aggregate moments\n",
    "\n",
    "Health inequalities: we use the gradient for the 4th quartile relative to the first. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4acc8877",
   "metadata": {},
   "outputs": [],
   "source": [
    "grad_c  = [1.061, 1.276]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0365aed2",
   "metadata": {},
   "source": [
    "GDP share of health expenditures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ac3df66c",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_c  = [0.0955, 0.1504]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b855728",
   "metadata": {},
   "source": [
    "Income"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "20c3dd9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_c  = [0.78, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa65a7cb",
   "metadata": {},
   "source": [
    "For health status, we use the average health from the ADL definition in Table 4. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "be08d6d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "h_c  = [0.931, 0.892]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a1eab219",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.096, 0.93 , 1.06 ],\n",
       "       [0.15 , 0.89 , 1.29 ]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi = np.zeros((Ncount,3))\n",
    "for ii in range(Ncount):\n",
    "    psi[ii,:] = [s_c[ii],h_c[ii],grad_c[ii]]\n",
    "\n",
    "psi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d01473d0",
   "metadata": {},
   "source": [
    "## Income inequalities: the US\n",
    "\n",
    "Table 5 gives the variance of the income process. If log income has mean 0 and variance equal to $\\sigma^2_{\\epsilon}$, we can solve for the quartiles (we use mid-points of quartiles, 12.5, 37.5, 62.5 and 87.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48a95d99",
   "metadata": {},
   "source": [
    "We take the mean for EU in Table 5. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e4046ef5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.19111728571428574, 0.486429]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sig_eu  = np.mean([0.256367, 0.094643, 0.214538, 0.237439, 0.169834, 0.094643, 0.270357])\n",
    "sig_c   = [sig_eu, 0.486429]\n",
    "sig_c"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e5a0143",
   "metadata": {},
   "source": [
    "We find the percentiles, transform in levels and normalize to 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a9f5396e",
   "metadata": {},
   "outputs": [],
   "source": [
    "qs   = [0.125,0.375,0.625,0.875]\n",
    "yy_c = np.zeros((Ncount,4))\n",
    "\n",
    "for ii in range(Ncount):\n",
    "    yy    = [norm(0,np.sqrt(sig_c[ii])).ppf(q) for q in qs]\n",
    "    yy    = np.exp(yy)\n",
    "    yy_c[ii,:] = y_c[ii]*yy/np.mean(yy)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9017749",
   "metadata": {},
   "source": [
    "Here are relative incomes. Much more inequality in US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "806c601a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4.975918497985015, 2.7340817164336597)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "yy_c[-1,-1]/yy_c[-1,0], yy_c[0,-1]/yy_c[0,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f2a9f4f",
   "metadata": {},
   "source": [
    "## Solution for the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "95c2f74f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/flangot/opt/anaconda3/lib/python3.7/site-packages/scipy/optimize/minpack.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last ten iterations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([7.73297119, 2.93571971, 1.43153705])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_c   = (psi[1,0], psi[1,1], psi[1,2], 1, yy_c[1], sigma)\n",
    "pd_c0     = [6 , 2 , 2 ] \n",
    "calibsol  = fsolve(hetero_calib, pd_c0, args=param_c) \n",
    "calibsol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8945fb87",
   "metadata": {},
   "source": [
    "We check the solution: to find the health gradient with the model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "214a0e9a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.        , 0.08376704, 0.18351354, 0.33271942]), 1.290000000028967)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv = np.zeros(yy.size)\n",
    "\n",
    "for ii in range(yy.size):\n",
    "    param_m   = (1, yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolv[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "msolv, ( 1-np.exp(-(calibsol[2]+calibsol[0]*msolv[-1])) )/( 1-np.exp(-(calibsol[2]+calibsol[0]*msolv[0])) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1543e6dd",
   "metadata": {},
   "source": [
    "## Solution for the Europe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "986c98f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calibsol_eu = np.zeros((1,2))\n",
    "param_s           = (psi[0,0], psi[0,1], yy_c[0], sigma, calibsol[1], calibsol[0])\n",
    "pd_c0             = [0.9 , calibsol[2]]\n",
    "calibsol_eu       = fsolve(hetero_calib_eu, pd_c0, args=param_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "50d82726",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.57669416, 1.86046548])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calibsol_eu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d3a89f8",
   "metadata": {},
   "source": [
    "## Individual decisions and health gradient: Europe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "dd23a4f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.01690788, 0.093706  , 0.15885468, 0.24990552]), 1.132026772925869)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv_eu = np.zeros(yy.size)\n",
    "for jj in range(yy.size):\n",
    "    param_m         = (calibsol_eu[0], yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eu[jj]    = max(0,fsolve(medexp, m0, args=param_m))\n",
    "grads_eu = ( 1-np.exp(-(calibsol_eu[1]+calibsol[0]*msolv_eu[-1])) )/(1-np.exp(-(calibsol_eu[1]+calibsol[0]*msolv_eu[0])))\n",
    "\n",
    "    \n",
    "msolv_eu, grads_eu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "597eefdf",
   "metadata": {},
   "source": [
    "## Demand : price elasticity in the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f3322eb8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.2792993863472035,\n",
       " 0.3384588537552355,\n",
       " 0.7350457180497365,\n",
       " 0.6551823760479963)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolvp = np.zeros(yy.size)\n",
    "msolvm = np.zeros(yy.size)\n",
    "\n",
    "for ii in range(yy.size):\n",
    "    param_m   = (1.1, yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvp[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m   = (0.9, yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvm[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "dpop    = (1.1-1)/1\n",
    "share0  = s_c[1]\n",
    "share1p = 1.1*np.mean(msolvp)\n",
    "share1m = 0.9*np.mean(msolvm)\n",
    "level0  = s_c[1]\n",
    "level1p = np.mean(msolvp)\n",
    "level1m = np.mean(msolvm)\n",
    "\n",
    "((share1p-share0)/share0)/dpop, -((share1m-share0)/share0 )/dpop, ((level1m-level0)/level0)/dpop, -((level1p-level0)/level0 )/dpop\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "1ac7c9c0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.30887912005121954, 0.6951140470488664)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((((share1p-share0)/share0)/dpop)+(-((share1m-share0)/share0)/dpop))/2, ((((level1m-level0)/level0)/dpop)+(-((level1p-level0)/level0)/dpop))/2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2db395ff",
   "metadata": {},
   "source": [
    "In the U.S., the price elasticity of the GDP share of health expenditures ($s_{US} = p_{US}m_{US}/y_{US}$) is equal to the price elasticity of the health expenditures ($p_{US}m_{US}$) becuase $y_{US}=1$.\n",
    "\n",
    "But, we can compute the price elasticity of the expenditures ($m_{US}$)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c4bbc61",
   "metadata": {},
   "source": [
    "## Demand : income elasticity in the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "92110032",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.13781097412077215,\n",
       " 0.26640725605422205,\n",
       " 1.1515920715328485,\n",
       " 1.2397665304487981)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolvp = np.zeros(yy.size)\n",
    "msolvm = np.zeros(yy.size)\n",
    "\n",
    "for ii in range(yy.size):\n",
    "    param_m   = (1, 1.1*yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvp[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m   = (1, 0.9*yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvm[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "dyoy    = (1.1-1)/1\n",
    "share0  = s_c[1]\n",
    "share1p = np.mean(msolvp)/1.1\n",
    "share1m = np.mean(msolvm)/0.9\n",
    "level0  = s_c[1]\n",
    "level1p = np.mean(msolvp)\n",
    "level1m = np.mean(msolvm)\n",
    "\n",
    "((share1p-share0)/share0)/dyoy, -((share1m-share0)/share0)/dyoy, ((level1p-level0)/level0)/dyoy, -((level1m-level0)/level0)/dyoy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "f964f2f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.2021091150874971, 1.1956793009908233)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((((share1p-share0)/share0)/dyoy)+(-((share1m-share0)/share0)/dyoy))/2, ((((level1p-level0)/level0)/dyoy)+(-((level1m-level0)/level0)/dyoy))/2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3af07c7a",
   "metadata": {},
   "source": [
    "## Demand : price elasticity in EU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ba949ad1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.07653415078077341,\n",
       " 0.0013215520589766755,\n",
       " 1.0260731657991429,\n",
       " 0.910292320053614)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv_eup = np.zeros(yy.size)\n",
    "msolv_eum = np.zeros(yy.size)\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m         = (calibsol_eu[0]*1.1, yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eup[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m         = (calibsol_eu[0]*0.9, yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eum[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "    \n",
    "share0  = s_c[0] \n",
    "share1p = calibsol_eu[0]*1.1*np.mean(msolv_eup)/np.mean(yy_c[0,:])\n",
    "share1m = calibsol_eu[0]*0.9*np.mean(msolv_eum)/np.mean(yy_c[0,:])\n",
    "\n",
    "level0  = s_c[0]/(calibsol_eu[0]/y_c[0]) \n",
    "level1p = np.mean(msolv_eup) \n",
    "level1m = np.mean(msolv_eum) \n",
    "\n",
    "((share1m-share0)/share0)/dpop, -((share1p-share0)/share0)/dpop, ((level1m-level0)/level0)/dpop, -((level1p-level0)/level0)/dpop \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8e41e4a9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.03760629936089837, 0.9681827429263785)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((((share1m-share0)/share0)/dpop)+(-((share1p-share0)/share0)/dpop))/2, ((((level1m-level0)/level0)/dpop)+(-((level1p-level0)/level0)/dpop))/2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdaae432",
   "metadata": {},
   "source": [
    "## Demand : income elasticity in EU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "ee03d8bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.6460609880115763, 0.864507865703855, 1.7106670868127356, 1.7780570791334676)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv_eup = np.zeros(yy.size)\n",
    "msolv_eum = np.zeros(yy.size)\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m         = (calibsol_eu[0], 1.1*yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eup[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m         = (calibsol_eu[0], 0.9*yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eum[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "share0  = s_c[0] \n",
    "share1p = calibsol_eu[0]*np.mean(msolv_eup)/(1.1*np.mean(yy_c[0,:]))\n",
    "share1m = calibsol_eu[0]*np.mean(msolv_eum)/(0.9*np.mean(yy_c[0,:]))\n",
    "    \n",
    "level0  = s_c[0]/(calibsol_eu[0]/y_c[0])    \n",
    "level1p = np.mean(msolv_eup) \n",
    "level1m = np.mean(msolv_eum)\n",
    "\n",
    "((share1p-share0)/share0)/dyoy, -((share1m-share0)/share0)/dyoy, ((level1p-level0)/level0)/dyoy, -((level1m-level0)/level0)/dyoy \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "1e6bf63d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.7552844268577157, 1.7443620829731015)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((((share1p-share0)/share0)/dpop)+(-((share1m-share0)/share0)/dpop))/2, ((((level1p-level0)/level0)/dpop)+(-((level1m-level0)/level0)/dpop))/2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf13ddbd",
   "metadata": {},
   "source": [
    "# Forecasting of the GDP share of health expenditures\n",
    "\n",
    "Observation: the shift in GDP per capita\n",
    "\n",
    "Implication: the implicit shift in healthcare prices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5fa8267",
   "metadata": {},
   "source": [
    "### GDP per capita\n",
    "\n",
    "GDP share of health expenditures: OECD data\n",
    "\n",
    "Real GDP per capita: WB data (World Development Indicators)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "2d32641c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1970\n",
    "GDPpc_eu70 = [17887.10176179620, 26927.95903692180, 17519.42717022000, 15729.85644157820, 21354.25810368360, 24388.29821013090, 11431.59455874400]\n",
    "PMoY_eu70  = [5.713, 7.676666666666667,   5.184, float(\"NAN\"), 5.434246479956248,  5.506, 3.145]\n",
    "GDPpc70    = [np.nanmean(GDPpc_eu70), 25279.18879113590]\n",
    "PMoY70     = [np.nanmean(PMoY_eu70) , 6.232]\n",
    "# 1990\n",
    "GDPpc_eu90 = [29473.85438511120,   39295.24344249190,   28512.10870941340,   27479.57862963770,   31093.54486501490,   34570.30573399390,   18962.44605133250]\n",
    "PMoY_eu90  = [8.031, 8.038, 7.958, 7.015, 7.082, 7.256, 6.101]\n",
    "GDPpc90    = [np.nanmean(GDPpc_eu90) , 39278.59534978750]\n",
    "PMoY90     = [np.nanmean(PMoY_eu90) , 11.273]\n",
    "\n",
    "# 2000\n",
    "GDPpc_eu00 = [34476.20802831070, 49241.92800413450, 33583.85744208620, 32337.89674255510, 40440.67534426340, 41117.15243035039, 23928.67165633810]\n",
    "PMoY_eu00  = [9.827999999999999, 8.103999999999999, 9.541, 7.58, 7.056, 7.412, 6.8160]\n",
    "GDPpc00 = [np.nanmean(GDPpc_eu00) , 48689.03128985530]\n",
    "PMoY00  = [np.nanmean(PMoY_eu00) , 12.502]\n",
    "\n",
    "# 2017\n",
    "GDPpc_eu10 = [42622.40993268510, 55735.76490114850, 37678.92729753470, 31231.66444646460, 46978.44880120770, 52576.81028210510, 27213.45166906160]\n",
    "PMoY_eu10  = [11.272, 10.218, 11.458, 8.9010, 10.143, 10.916, 8.840999999999999]\n",
    "GDPpc10    = [np.nanmean(GDPpc_eu10) , 58387.77580835719]\n",
    "PMoY10     = [np.nanmean(PMoY_eu10) , 17.149999999999999]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "4e128ceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "psi0    = [i/100 for i in PMoY90]\n",
    "psi1    = [i/100 for i in PMoY10]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba04b0a3",
   "metadata": {},
   "source": [
    "## Income heterogeneity by country (data)\n",
    "We rescale the income level and compute the quartiles in 1970 and 2017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "aab01b00",
   "metadata": {},
   "outputs": [],
   "source": [
    "y70_c = y_c*np.array(GDPpc70)/np.array(GDPpc00)\n",
    "y90_c = y_c*np.array(GDPpc90)/np.array(GDPpc00)\n",
    "y00_c = y_c*np.array(GDPpc00)/np.array(GDPpc00)\n",
    "y10_c = y_c*np.array(GDPpc10)/np.array(GDPpc00)\n",
    "\n",
    "yy70_c = np.zeros((Ncount,yy.size))\n",
    "yy90_c = np.zeros((Ncount,yy.size))\n",
    "yy00_c = np.zeros((Ncount,yy.size))\n",
    "yy10_c = np.zeros((Ncount,yy.size))\n",
    "\n",
    "for ii in range(Ncount):\n",
    "    yy    = [norm(0,np.sqrt(sig_c[ii])).ppf(q) for q in qs]\n",
    "    yy    = np.exp(yy)\n",
    "    yy70_c[ii,:] = y70_c[ii]*yy/np.mean(yy)\n",
    "    yy90_c[ii,:] = y90_c[ii]*yy/np.mean(yy)\n",
    "    yy00_c[ii,:] = y00_c[ii]*yy/np.mean(yy)\n",
    "    yy10_c[ii,:] = y10_c[ii]*yy/np.mean(yy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "de283a97",
   "metadata": {},
   "outputs": [],
   "source": [
    "prevsol0 = [0,0]\n",
    "prevsol1 = [0,0]\n",
    "\n",
    "#                  pm/y   income     sigma  phi          alpha_1         alpha_0\n",
    "param_s        = (psi0[0], yy90_c[0], sigma, calibsol[1], calibsol_eu[1], calibsol[0])\n",
    "pd_c0          = .1 #.3\n",
    "prevsol0[0]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n",
    "param_s        = (psi0[1], yy90_c[1], sigma, calibsol[1], calibsol[2], calibsol[0])\n",
    "pd_c0          = .8\n",
    "prevsol0[1]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n",
    "\n",
    "param_s        = (psi1[0], yy10_c[0], sigma, calibsol[1], calibsol_eu[1], calibsol[0])\n",
    "pd_c0          = .5#.5\n",
    "prevsol1[0]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n",
    "param_s        = (psi1[1], yy10_c[1], sigma, calibsol[1], calibsol[2], calibsol[0])\n",
    "pd_c0          = 1.2\n",
    "prevsol1[1]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "43bea0f2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([[0.2397345922351133], [0.48417706539364735]],\n",
       " [[0.5253143630056939], [1.3889148216652534]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prevsol0, prevsol1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "825acd69",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAATJklEQVR4nO3dYYhdZ37f8e9vZYtdSKnsemyEJCqzDKFqIFohtIJ9s6yzRdKWjvfFtlaJLYxBFrFgAyllmr6ok1fONhsXg5CQuyIyXSIMSfHgTCOMuiYsVI7GW6/WWkd4Itz1rAd74rBOjEOM7H9ezFFze3ule0Yaa9Z6vh843HOe5//c+xy43J/Oo3PnpqqQJLXnM2s9AUnS2jAAJKlRBoAkNcoAkKRGGQCS1Kjb1noCK3HXXXfV1q1b13oakvSp8vLLL/9lVU0Mt3+qAmDr1q3Mzc2t9TQk6VMlyf8Z1e4SkCQ1ygCQpEYZAJLUKANAkhplAEhSo3oFQJI9SS4mmU8yPaI/SZ7q+s8n2THUvy7J/07y/EDbnUleSPJ693jHjZ+OJKmvsQGQZB1wBNgLbAP2J9k2VLYXmOy2g8DRof5vAq8NtU0DZ6pqEjjTHUuSbpI+VwC7gPmqulRVHwKngKmhmingmVp2FtiQZCNAks3A14D/OmLMyW7/JHD/9Z2CJOl69AmATcCbA8cLXVvfmv8C/Hvg46Ex91TVIkD3ePeoF09yMMlckrmlpaUe05Uk9dHnm8AZ0Tb8KzIja5L8S+Cdqno5yZdXOLflJ6k6DhwH2Llzp79eo1vW1uk/Xusp6OfYG098bdWfs88VwAKwZeB4M/BWz5ovAf8qyRssLx19Jcl/62reHlgm2gi8s+LZS5KuW58AOAdMJrk3yXrgAWBmqGYGeKi7G2g38F5VLVbVf6iqzVW1tRv3P6vqVwfGHOj2DwDP3ejJSJL6G7sEVFWXkxwGTgPrgBNVdSHJoa7/GDAL7APmgQ+Ah3u89hPAs0keAX4CfOP6TkGSdD16/TXQqppl+UN+sO3YwH4Bj415jheBFweO3wXu6z9VSdJq8pvAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1KheAZBkT5KLSeaTTI/oT5Knuv7zSXZ07Z9N8mdJfpjkQpLfGhjzeJKfJnml2/at3mlJksYZ+5OQSdYBR4CvAgvAuSQzVfXjgbK9wGS3fRE42j3+HfCVqno/ye3A95P8j6o62417sqp+d/VOR5LUV58rgF3AfFVdqqoPgVPA1FDNFPBMLTsLbEiysTt+v6u5vdtqtSYvSbp+fQJgE/DmwPFC19arJsm6JK8A7wAvVNVLA3WHuyWjE0nuGPXiSQ4mmUsyt7S01GO6kqQ++gRARrQN/yv+qjVV9VFVbQc2A7uS/FLXfxT4PLAdWAS+PerFq+p4Ve2sqp0TExM9pitJ6qNPACwAWwaONwNvrbSmqn4GvAjs6Y7f7sLhY+BplpeaJEk3SZ8AOAdMJrk3yXrgAWBmqGYGeKi7G2g38F5VLSaZSLIBIMnngF8B/rw73jgw/uvAqzd2KpKklRh7F1BVXU5yGDgNrANOVNWFJIe6/mPALLAPmAc+AB7uhm8ETnZ3En0GeLaqnu/6vpVkO8tLRW8Aj67WSUmSxhsbAABVNcvyh/xg27GB/QIeGzHuPPCFqzzngyuaqSRpVflNYElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWpUrwBIsifJxSTzSaZH9CfJU13/+SQ7uvbPJvmzJD9MciHJbw2MuTPJC0le7x7vWL3TkiSNMzYAut/zPQLsBbYB+5NsGyrbC0x220HgaNf+d8BXquqXge3Anu5H4wGmgTNVNQmc6Y4lSTdJnyuAXcB8VV2qqg+BU8DUUM0U8EwtOwtsSLKxO36/q7m922pgzMlu/yRw/w2chyRphfoEwCbgzYHjha6tV02SdUleAd4BXqiql7qae6pqEaB7vHvUiyc5mGQuydzS0lKP6UqS+ugTABnRVn1rquqjqtoObAZ2JfmllUywqo5X1c6q2jkxMbGSoZKka7itR80CsGXgeDPw1kprqupnSV4E9gCvAm93y0SLSTayfIXwidk6/cef5NPrU+yNJ7621lOQ1kSfK4BzwGSSe5OsBx4AZoZqZoCHuruBdgPvdR/sE0k2ACT5HPArwJ8PjDnQ7R8AnruxU5EkrcTYK4CqupzkMHAaWAecqKoLSQ51/ceAWWAfMA98ADzcDd8InOzuJPoM8GxVPd/1PQE8m+QR4CfAN1bvtCRJ4/RZAqKqZln+kB9sOzawX8BjI8adB75wled8F7hvJZOVJK0evwksSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjeoVAEn2JLmYZD7J9Ij+JHmq6z+fZEfXviXJ95K8luRCkm8OjHk8yU+TvNJt+1bvtCRJ44z9Scju93yPAF8FFoBzSWaq6scDZXuByW77InC0e7wM/EZV/SDJPwJeTvLCwNgnq+p3V+90JEl99bkC2AXMV9WlqvoQOAVMDdVMAc/UsrPAhiQbq2qxqn4AUFV/A7wGbFrF+UuSrlOfANgEvDlwvMD//yE+tibJVpZ/IP6lgebD3ZLRiSR3jHrxJAeTzCWZW1pa6jFdSVIffQIgI9pqJTVJfgH4Q+DXq+qvu+ajwOeB7cAi8O1RL15Vx6tqZ1XtnJiY6DFdSVIffQJgAdgycLwZeKtvTZLbWf7w/25V/dGVgqp6u6o+qqqPgadZXmqSJN0kfQLgHDCZ5N4k64EHgJmhmhngoe5uoN3Ae1W1mCTAd4DXqur3Bgck2Thw+HXg1es+C0nSio29C6iqLic5DJwG1gEnqupCkkNd/zFgFtgHzAMfAA93w78EPAj8KMkrXdtvVtUs8K0k21leKnoDeHSVzkmS1MPYAADoPrBnh9qODewX8NiIcd9n9P8PUFUPrmimkqRV5TeBJalRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVG9AiDJniQXk8wnmR7RnyRPdf3nk+zo2rck+V6S15JcSPLNgTF3Jnkhyevd4x2rd1qSpHHGBkCSdcARYC+wDdifZNtQ2V5gstsOAke79svAb1TVPwN2A48NjJ0GzlTVJHCmO5Yk3SR9rgB2AfNVdamqPgROAVNDNVPAM7XsLLAhycaqWqyqHwBU1d8ArwGbBsac7PZPAvff2KlIklaiTwBsAt4cOF7gHz7Ee9ck2Qp8AXipa7qnqhYBuse7R714koNJ5pLMLS0t9ZiuJKmPPgGQEW21kpokvwD8IfDrVfXX/acHVXW8qnZW1c6JiYmVDJUkXUOfAFgAtgwcbwbe6luT5HaWP/y/W1V/NFDzdpKNXc1G4J2VTV2SdCP6BMA5YDLJvUnWAw8AM0M1M8BD3d1Au4H3qmoxSYDvAK9V1e+NGHOg2z8APHfdZyFJWrHbxhVU1eUkh4HTwDrgRFVdSHKo6z8GzAL7gHngA+DhbviXgAeBHyV5pWv7zaqaBZ4Ank3yCPAT4BurdlaSpLHGBgBA94E9O9R2bGC/gMdGjPs+o/9/gKp6F7hvJZOVJK0evwksSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjeoVAEn2JLmYZD7J9Ij+JHmq6z+fZMdA34kk7yR5dWjM40l+muSVbtt346cjSeprbAAkWQccAfYC24D9SbYNle0FJrvtIHB0oO/3gT1Xefonq2p7t81epUaS9AnocwWwC5ivqktV9SFwCpgaqpkCnqllZ4ENSTYCVNWfAn+1mpOWJN24PgGwCXhz4Hiha1tpzSiHuyWjE0nuGFWQ5GCSuSRzS0tLPZ5SktRHnwDIiLa6jpphR4HPA9uBReDbo4qq6nhV7ayqnRMTE2OeUpLUV58AWAC2DBxvBt66jpr/R1W9XVUfVdXHwNMsLzVJkm6SPgFwDphMcm+S9cADwMxQzQzwUHc30G7gvapavNaTXvk/gs7XgVevVitJWn23jSuoqstJDgOngXXAiaq6kORQ138MmAX2AfPAB8DDV8Yn+QPgy8BdSRaA/1RV3wG+lWQ7y0tFbwCPrt5pSZLGGRsAAN0tmrNDbccG9gt47Cpj91+l/cH+05QkrTa/CSxJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmN6hUASfYkuZhkPsn0iP4kearrP59kx0DfiSTvJHl1aMydSV5I8nr3eMeNn44kqa+xAZBkHXAE2AtsA/Yn2TZUtheY7LaDwNGBvt8H9ox46mngTFVNAme6Y0nSTdLnCmAXMF9Vl6rqQ+AUMDVUMwU8U8vOAhuSbASoqj8F/mrE804BJ7v9k8D91zF/SdJ16hMAm4A3B44XuraV1gy7p6oWAbrHu0cVJTmYZC7J3NLSUo/pSpL66BMAGdFW11FzXarqeFXtrKqdExMTq/GUkiT6BcACsGXgeDPw1nXUDHv7yjJR9/hOj7lIklZJnwA4B0wmuTfJeuABYGaoZgZ4qLsbaDfw3pXlnWuYAQ50+weA51Ywb0nSDRobAFV1GTgMnAZeA56tqgtJDiU51JXNApeAeeBp4NeujE/yB8D/An4xyUKSR7quJ4CvJnkd+Gp3LEm6SW7rU1RVsyx/yA+2HRvYL+Cxq4zdf5X2d4H7es9UkrSq/CawJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNapXACTZk+Rikvkk0yP6k+Sprv98kh3jxiZ5PMlPk7zSbftW55QkSX2MDYAk64AjwF5gG7A/ybahsr3AZLcdBI72HPtkVW3vtlkkSTdNnyuAXcB8VV2qqg+BU8DUUM0U8EwtOwtsSLKx51hJ0hroEwCbgDcHjhe6tj4148Ye7paMTiS5Y9SLJzmYZC7J3NLSUo/pSpL66BMAGdFWPWuuNfYo8HlgO7AIfHvUi1fV8araWVU7JyYmekxXktTHbT1qFoAtA8ebgbd61qy/2tiqevtKY5Knged7z1qSdMP6XAGcAyaT3JtkPfAAMDNUMwM81N0NtBt4r6oWrzW2+z+CK74OvHqD5yJJWoGxVwBVdTnJYeA0sA44UVUXkhzq+o8Bs8A+YB74AHj4WmO7p/5Wku0sLwm9ATy6iuclSRqjzxIQ3S2as0Ntxwb2C3is79iu/cEVzVSStKr8JrAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1qlcAJNmT5GKS+STTI/qT5Kmu/3ySHePGJrkzyQtJXu8e71idU5Ik9TE2AJKsA44Ae4FtwP4k24bK9gKT3XYQONpj7DRwpqomgTPdsSTpJulzBbALmK+qS1X1IXAKmBqqmQKeqWVngQ1JNo4ZOwWc7PZPAvff2KlIklaiz4/CbwLeHDheAL7Yo2bTmLH3VNUiQFUtJrl71IsnOcjyVQXA+0ku9pizxrsL+Mu1nsTPg/zOWs9AV+F7dMANvk//6ajGPgGQEW3Vs6bP2GuqquPA8ZWM0XhJ5qpq51rPQ7oa36OfvD5LQAvAloHjzcBbPWuuNfbtbpmI7vGd/tOWJN2oPgFwDphMcm+S9cADwMxQzQzwUHc30G7gvW5551pjZ4AD3f4B4LkbPBdJ0gqMXQKqqstJDgOngXXAiaq6kORQ138MmAX2AfPAB8DD1xrbPfUTwLNJHgF+AnxjVc9M47ispp93vkc/Yala0ZK8JOkW4TeBJalRBoAkNcoAuMUl+SjJKwPbdNf+RpK7Buq+nOT5tZupWpVka5JXh9oeT/LvkuxO8lL33n0tyeNrNM1bUp/vAejT7W+ravtaT0K6TieBf11VP+z+tMwvrvWEbiUGgKSfZ3cDV/5iwEfAj9d2OrcWl4BufZ8bWgL6N2s9IWkFngQuJvnvSR5N8tm1ntCtxCuAW9/VloBG3f/rPcFaC1d731VV/XaS7wL/Avi3wH7gyzdrYrc6rwDa9S4w+BsMd+If3tLaGH4vwsD7sar+oqqOAvcBv5zkn9zk+d2yDIB2vQg8CP/3dxt+FfjeWk5Ibaqq94HFJPfB8o9FAXuA7yf5WpIrf1RyEvgI+NmaTPQW5DeBb3FJPgJ+NND0J1U1neQfs/zDPf+c5b/a+ifAdFV9vAbTVOO6H4o6wj9cCfznqvpuklPADpb/xMxl4D9W1ek1muYtxwCQpEa5BCRJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqP+Htx6+zX1X9E6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "columns = ('EU', 'US')\n",
    "Price0 = sum(prevsol0,[])\n",
    "Price1 = sum(prevsol1,[]) \n",
    "gPrice = [((Price1[ii]/Price0[ii])**(1/(2017-1990)))-1 for ii in np.arange(len(columns))]\n",
    "plt.bar(np.arange(len(columns)), gPrice)\n",
    "plt.xticks(np.arange(len(columns)), labels=columns)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "9bd05d30",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.029480429856935553, 0.0398023450155367]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gPrice"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23335f30",
   "metadata": {},
   "source": [
    "### Size of the error when price changes are omitted\n",
    "\n",
    "If prices change, then the model predicts exactly the GDP share of health expenditures (targeted data)\n",
    "\n",
    "$\\Rightarrow$ If prices do not change, then it exists a gap between model and data: this gap gives the error if price change are neglected "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "a119d338",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([8.165492112286433, 13.969329123513006],\n",
       " [7.3544285714285715, 11.273],\n",
       " [10.463701949155052, 15.478811642261823],\n",
       " [10.249857142857142, 17.15])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv0_90 = np.zeros((2,yy.size))\n",
    "pmoys0_90 = [0,0]\n",
    "msolv0_10 = np.zeros((2,yy.size))\n",
    "pmoys0_10 = [0,0]\n",
    "\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m        = (calibsol_eu[0], yy90_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0             = .1\n",
    "    msolv0_90[0,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m        = (calibsol_eu[0], yy10_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0             = .1\n",
    "    msolv0_10[0,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m        = (1, yy90_c[1,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0             = .1\n",
    "    msolv0_90[1,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m        = (1, yy10_c[1,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0             = .1\n",
    "    msolv0_10[1,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "price_bench = [calibsol_eu[0],1] \n",
    "for ii in range(2):            \n",
    "    pmoys0_90[ii]= 100*price_bench[ii]*np.mean(msolv0_90[ii,:])/np.mean(yy90_c[ii,:])\n",
    "    pmoys0_10[ii]= 100*price_bench[ii]*np.mean(msolv0_10[ii,:])/np.mean(yy10_c[ii,:])\n",
    "\n",
    "pmoys0_90, PMoY90, pmoys0_10, PMoY10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83f75020",
   "metadata": {},
   "source": [
    "Without price change, the model over-estimate the GDP share of health expenditures in the 1990s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "2cdc6ddb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.8110635408578615, 2.6963291235130065]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[pmoys0_90[ii]-PMoY90[ii] for ii in np.arange(2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "267e1682",
   "metadata": {},
   "source": [
    "Without price changes, the model under-estimate the GDP share of health expenditures in the 2010s for the U.S."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "ef97b436",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.21384480629791014, -1.6711883577381759]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[pmoys0_10[ii]-PMoY10[ii] for ii in np.arange(2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef478f4c",
   "metadata": {},
   "source": [
    "# Forcasting long run change in GDP share of Health Expenditures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7df75ef3",
   "metadata": {},
   "source": [
    "### Long run changes in GDP per capita and income dispersion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "cc424891",
   "metadata": {},
   "outputs": [],
   "source": [
    "gGDPpc    = [((GDPpc10[ii]/GDPpc90[ii])**(1/(2017-1990))-1) for ii in np.arange(2)]\n",
    "agGDPpc   = np.array(gGDPpc)\n",
    "agGDPpc10 = np.array(GDPpc10)\n",
    "\n",
    "horizon = 4\n",
    "fGDPpc  = np.zeros((2,horizon))\n",
    "fy_c    = np.zeros((2,horizon))\n",
    "fyy_c   = np.zeros((2,horizon,yy.size))\n",
    "\n",
    "for jj in range(2):\n",
    "    for ii in range(horizon):\n",
    "        fGDPpc[jj,ii] = (GDPpc10[jj]*(1+gGDPpc[jj])**(2030+(ii*10)-2017))/GDPpc00[jj]\n",
    "        fy_c[jj,ii] = y_c[jj]*fGDPpc[jj,ii]\n",
    "        yy    = [norm(0,np.sqrt(sig_c[jj])).ppf(q) for q in qs]\n",
    "        yy    = np.exp(yy)\n",
    "        fyy_c[jj,ii,:] = fy_c[jj,ii]*yy/np.mean(yy)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "2d40c0d9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.35719765, 1.53905893, 1.74528919, 1.97915382],\n",
       "       [1.45139469, 1.68093419, 1.94677557, 2.25466004]])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fGDPpc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "2b8441b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.01265430880932894, 0.014790792069827585]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gGDPpc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a307dd87",
   "metadata": {},
   "source": [
    "### Long run changes in prices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "0f77bc8f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.76639827, 1.02479138, 1.3703024 , 1.83230332],\n",
       "       [2.30693796, 3.40834729, 5.03560628, 7.43977312]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fPrice  = np.zeros((2,horizon))\n",
    "\n",
    "for jj in range(2):\n",
    "    for ii in range(horizon):\n",
    "        fPrice[jj,ii] = (Price1[jj]*(1+gPrice[jj])**(2030+(ii*10)-2017))\n",
    "        \n",
    "fPrice"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e29d255",
   "metadata": {},
   "source": [
    "Europe with changes in price"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "bbfbdab0",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx     = np.zeros((horizon,yy.size))\n",
    "yy_eu  = np.zeros(horizon)\n",
    "yy_eu0 = np.zeros(horizon)\n",
    "yy_us  = np.zeros(horizon)\n",
    "yy_us0 = np.zeros(horizon)\n",
    "\n",
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (fPrice[0,ii], fyy_c[0,ii,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "        m0               = .05#.01\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_eu[ii] = 100*fPrice[0,ii]*np.mean(xx[ii,:])/np.mean(fyy_c[0,ii,:])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c852a7b4",
   "metadata": {},
   "source": [
    "Europe without price changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "ce3c0840",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (Price1[0], fyy_c[0,ii,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "        m0               = .01#.005\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_eu0[ii] = 100*Price1[0]*np.mean(xx[ii,:])/np.mean(fyy_c[0,ii,:])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99c5497f",
   "metadata": {},
   "source": [
    "The U.S. with price changes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ce7da493",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (fPrice[1,ii], fyy_c[1,ii,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "        m0               = .01#.005\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_us[ii] = 100*fPrice[1,ii]*np.mean(xx[ii,:])/np.mean(fyy_c[1,ii,:])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2079c7d",
   "metadata": {},
   "source": [
    "The U.S. without price changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "58fb5606",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = np.zeros((horizon,yy.size))\n",
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (Price1[1], fyy_c[1,ii,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "        m0               = .01#.005\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_us0[ii] = 100*Price1[1]*np.mean(xx[ii,:])/np.mean(fyy_c[1,ii,:])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "8003a56b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD4CAYAAAD2FnFTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8kUlEQVR4nO3de3xU9Zn48c8394SEAAkCSYBAEgJIFRVETYAIEQGRoLagrRfUitpibeuv3Wrbl65bd1lr27VLQRGs4gpKBQKyBEVkFfAKFOUSQghECQmXcAuBhCQzz++PmYyTZCaZJDOZXJ7365XXzJw558zDl+Q88z3ne76PERGUUkqp+gL8HYBSSqn2SROEUkoplzRBKKWUckkThFJKKZc0QSillHIpyN8BeFNsbKwkJib6OwyllOowduzYUSoivV2916kSRGJiItu3b/d3GEop1WEYY75x956eYlJKKeWSJgillFIuaYJQSinlUqe6BuFKdXU1RUVFVFZW+juUTiksLIyEhASCg4P9HYpSyss6fYIoKioiKiqKxMREjDH+DqdTERFOnTpFUVERgwYN8nc4Sikv6/QJorKyUpODjxhjiImJ4eTJk/4ORakuJy0tjdLS0gbLY2Nj2bZtm1c+o0tcg9Dk4Dvatkr5h6vk0NjylugSCUIppVTz+ewUkzGmP7AU6AtYgUUi8qIx5m0g1b5aD+CsiIx0sX0hcB6wADUiMspXsSqlVEdRUVHBxx9/3Caf5ctrEDXAEyKy0xgTBewwxmwUkVm1Kxhj/gSca2QfN4qI9/pLSinVAVVWVrJlyxZycnLYvHkzFy9ebJPP9dkpJhEpEZGd9ufngVwgvvZ9Yzt5PRNY7qsY2pPAwEBGjhzp+Jk3bx6FhYWMGDGiznrPPPMML7zwQpP7W716NcYY9u/f76uQ2bBhA6mpqSQnJzNv3jyffY5SqqFLly7xwQcf8MQTT3D99dczd+5cPvnkE2699VZee+21NomhTUYxGWMSgauAz50WjwWOi0i+m80EeN8YI8DLIrLIzb7nAHMABgwY0OpYt27dyooVKygtLSU2NpaZM2eSnp7e6v2Gh4eza9euOssKCwtbvL/ly5czatQo3nrrLZ555plWxWaxWCguLqZ///51lv30pz9l48aNJCQkMHr0aKZPn87w4cNb9VlKKfeqqqocPYUPP/yQCxcu0KNHD6ZNm8aUKVO49tprCQqyHbZjY2PdjmLyFp8nCGNMJLAS+LmIlDm9dReN9x7SRKTYGHMZsNEYs19EGpx4syeORQCjRo1qVYHtrVu3snjxYqqqqgDbaIDFixcDeCVJeEt5eTkfffQRGzdu5Ac/+EGdBHHnnXciIhQWFnLs2DEWLFjALbfc0uj+nnjiCZKSknjsscccy7744guSk5MZPHiwY79r1qzRBKGUl1VVVbFt2zZycnLYtGkT5eXl9OjRgylTpjBlyhTGjBnj8kZUbw1lbYxPE4QxJhhbcnhTRFY5LQ8CbgeucbetiBTbH08YY1YD1wKtujKzdOlSvvnG7cSF5OfnU1NTU2dZVVUVixYtYvPmzS63GThwIPfee2+Tn11RUcHIkSMdr5988knGjBnjWeD1ZGdnk5mZyRVXXEG3bt3YuXMnV199NQBfffUVM2bM4O2332br1q388pe/bJAgMjMzOXbsGGC72S03N5fhw4czYMAAsrKyADh69GidHkVCQgKff/45SqnWq6qq4tNPPyUnJ4cPPviA8+fPEx0dzc0338yUKVO47rrr2sXsBL4cxWSAJUCuiPy53tuZwH4RKXKzbTcgQETO259PAp71Vay16ieHppY3h6tTTO6SVVP3Fixfvpw5c+YAMHPmTJYvX87VV19NRUUFpaWlPP300wAMHz6cM2fONNj+gw8+AOCTTz7h17/+NRcvXiQsLKzOOiINO2N6z4NSLVddXV0nKZSVlREVFcVNN93E5MmTuf766wkJCfF3mHX4sgeRBtwD7DbG7LIve0pE1gN3Uu/0kjEmDlgsIlOBPsBq+wEpCFgmIhtaG1BT3/R/9rOfuT2n9/vf/761H99ATExMgwP46dOnG5224tSpU3zxxResWmXrkM2aNYvx48fz/PPPs2fPHlJSUhwH+507d3LllVc22EdtD+LIkSN0796dUaNsI4ife+45Rw8iISGBI0eOOLYpKioiLi6udf9gpbqY6upqPvvsM0dSOHfuHJGRkWRmZjJlyhRuuOGGdpcUnPksQYjIVsDlV04Rme1iWTEw1f78ENDwyOZjM2fOrHMNAiAkJISZM2f65PMiIyPp168fmzZtYuLEiZw+fZoNGzbw+OOPAzBx4kSWLl1KfLxj8BfvvPMOU6dOJTQ0FIBBgwbRt29ftm7dSl5eHt9++y2VlZVYLBaefvppnn/++Qaf++6773LTTTexZMkSt6e5Ro8eTX5+PocPHyY+Pp633nqLZcuW+aAVlOpcampq+Pzzz8nJyWHjxo2cPXuWbt26MXHiRKZMmUJ6enq7TgrOOv1cTM1ReyF6xYoVnDp1ipiYGK+NYqp/DWLy5MnMmzePpUuX8tOf/pQnnngCgKeffpqkpCSsVisHDx6kV69edfazfPlyvv76a5xLq546dYply5YRFBTEj370IzIyMigrK+Opp54iLS3NZTzPPfdco9dAgoKCmD9/PjfffDMWi4UHHniAyy+/vOUNoFQnVlNTw5dffklOTg7vv/8+Z86cISIiok5SqP1S15EYV+eaO6pRo0ZJ/ZKjubm5DBs2zE8RtdyePXt49dVX+fOf61++cW/cuHG88sorpKamNr2yF3XUNlaqNSwWS52kcPr0aSIiIpgwYYIjKdS/ttceGWN2uJupQnsQ7dSIESOalRwACgoKSElJ8VFESimLxcKOHTvIycnhvffe49SpU4SHh3PjjTcyZcoUxo0b1yGSgqc0QXQiR48e9XcISnU6FouFnTt3OnoKJ0+eJCwsjIyMDKZMmcL48eMJDw/3d5g+oQlCKaXqsVqt7Ny5kw0bNrBhwwZOnjxJaGgoGRkZTJ48mYyMDCIiIvwdps9pglBKKWxJYdeuXeTk5LBhwwZOnDhBaGgo48ePd/QUunXr5u8w25QmCKVUlyUifPXVV46kcOzYMUJCQhg3bhxTpkwhIyODyMhIf4fpN5oglFJdiojw9ddfO04fFRcXExwczNixY3niiSeYMGFCl04KzjRBKKU6PRFh9+7djtFHR48eJTg4mPT0dB5//HEmTpxIVFSUv8NsdzRBKKU6JRFh79695OTkkJOT40gKaWlpPPbYY0ycOJHu3bv7O8x2TROEUqrTqJ2duDYpHDlyhKCgIG644Qbmzp3LxIkTiY6O9neYHYYmCKVUhyYi7N+/nw0bNpCTk8M333xDYGAg119/PY888giZmZn06NHD32F2SD4rOaq+01Rp0eeee47LL7+cK664gpEjR3pcd8HXZUe15Khqr2qTwl/+8hcmT57MjBkzeOWVV0hISOAPf/gDW7duZcmSJXz/+9/X5NAK2oNwkpaW5na6b19Vb/r0009Zt24dO3fuJDQ0lNLS0jqzyTbGW2VHteSo6ghEhPz8fMfpo8OHDxMQEMB1113HAw88wE033dRgckvVOpognLhKDo0t94aSkhJiY2MdMz16Wk/WXdlRLTmqOip3X9B69OjB3XffTU5ODgUFBQQEBHDttdcye/ZsbrrpJmJiYvwQbdfQpRLEc8891+LTMffcc4/L5UOHDuW3v/1ti2OaNGkSzz77LEOGDCEzM9NRAKgp7sqOaslR1VG5+yJ29uxZ/va3vzF69GjuvvtuJk2a5PEXKdU6XSpB+Iu7Up3GGCIjI9mxYwdbtmxh8+bNzJo1i3nz5jF79uxG9+mq7OiwYcO05KjqlLZs2ULv3r39HUaX06USRFPf9Buro/DGG2+0+HObKi0aGBhIRkYGGRkZfO973+P1119vNEG4Kzs6c+ZMLTmqOpTKyko+/PBDVq9e3eh6mhz8Q0cxtQHn0qKAo7Roeno6eXl55OfnO9bdtWsXAwcOBGwlR11N4e2u7OiGDRscJUcvXLjA008/zS9+8YsG27/77rv06NGD999/nyNHjrBnzx727NnjSA5Qt+RoVVUVb731FtOnT/dqu6iuSUTYvn07v/vd70hLS+MXv/gFBw4c8HdYygWf9SCMMf2BpUBfwAosEpEXjTHPAA8BJ+2rPiUi611sPxl4EQgEFouIz8dZxsbGuh3F1FruSovu2LGDxx57jLNnzxIUFERycjKLFi1yW3IU3Jcdveaaa7TkqGq3vv32W7Kzs1mzZg1FRUVERERw8803k5WVxZgxY7QqYTvks5Kjxph+QD8R2WmMiQJ2ADOAmUC5iLzQyLaBwAHgJqAI+BK4S0T2NfaZWnJUS46q9qWsrIycnByys7PZuXMnxhiuv/56srKymDRpUp2aCv4YZq78VHJUREqAEvvz88aYXCDew82vBQ6KyCEAY8xbQBbQaILoTLTkqOqoqqur2bZtG9nZ2WzatImqqiqSkpJ44oknmD59On379nW5nSaB9qdNLlIbYxKBq4DPgTRgrjHmXmA78ISI1B9qEw8ccXpdBLg8H2KMmQPMARgwYIB3A+9gtOSo8pfa4dLZ2dmsW7eOU6dO0bNnT2bNmsWMGTO4/PLLdRRcB+TzBGGMiQRWAj8XkTJjzELg3wCxP/4JeKD+Zi525fJcmIgsAhaB7RSTt+JWSjXt+PHjrFu3juzsbA4cOEBwcDATJkwgKyuLsWPHEhIS4u8QVSv4NEEYY4KxJYc3RWQVgIgcd3r/FWCdi02LgP5OrxOAYh+GqpTyUEVFBR988AHZ2dl88sknWK1WRo4cyTPPPMOUKVN07qNOxJejmAywBMgVkT87Le9nvz4BcBuwx8XmXwIpxphBwFHgTuCHvopVKdU4q9XK9u3bWb16Ne+99x4XLlwgLi6Ohx9+mKysLMc9Papz8WUPIg24B9htjNllX/YUcJcxZiS2U0aFwMMAxpg4bMNZp4pIjTFmLvAetmGur4rIXh/GqpRy4fDhw6xZs4a1a9dy9OhRIiIiHLOnjh49moAAvZWqM/PlKKatuL6W0OCeB/v6xcBUp9fr3a2rlPKds2fPsn79etasWcOuXbsICAjghhtu4Be/+AWZmZmEh4f7O0TVRrrUVBtKKdeqqqrYsmULa9as4cMPP6S6upqUlBR+9atfceutt9KnTx9/h6j8QBOEUl2UiLBnzx7WrFnDunXrOHPmDL169eKHP/whM2bMYNiwYTo0tYvTBKFUF3Ps2DHWrl3LmjVrOHjwICEhIUyYMIHbbruNtLQ0goOD/R2iaif0ClMbCQwMZOTIkY6fefPmNVmKtDG+LjcKWnK0M7lw4QLZ2dncf//9ZGRk8Kc//Ynu3bvz7LPPsnXrVl588UUyMjI0Oag6tAfhwqVLl9i3bx/Dhw93zJjaWuHh4ezatavOssLCwhbvz1vlRt3RkqMdn9Vq5fPPPyc7O5v333+fixcvkpCQwE9/+lOmT5/umDVYKXc0QbhQWFjIuXPnKCwsbPOJ7zzhrtwotKzkqCtacrTjKigocAxNLSkpITIykltuuYUZM2ZwzTXX6HUF5bEulSDy8/MpLy93+/65c+fqvC4pKaGkxHZPX3R0tMttIiMjPZogr6KigpEjRzpeP/nkk41Ot90Yd+VGAY9Kjo4dO5bz58832O8LL7xAZmYmoCVHO5ozZ86wfv16srOz+frrrwkICCA9PZ1f//rXTJgwoUHFQKU80aUSRFOioqKorKykurrasSw4ONgrf1yuTjF98803Ltdt6hueq3KjV199NRUVFR6VHN2yZUuT8WrJ0favqqqKjz76iOzsbD766COqq6sZOnQov/nNb5g2bZpWYVOt1qUShCff9PPy8igpKSEgIACr1UpsbKzPTjM1VYrUFXflRp9//nn27NnjUclRT3oQWnK0fRIRdu/eTXZ2Nv/7v//L2bNniY2N5Z577iErK4uhQ4f6O0TViXSpBOGJ6upq4uLiiIuLo7i4mKqqKp99lnMp0okTJzpKkT7++OOAreTo0qVLiY//royGu3KjW7duJS8vz1Fy1GKx8PTTT/P88883+FxPehDOJUfj4+N56623WLZsmZf+5aq5iouLWbt2LdnZ2Rw+fJjQ0FAyMzPJysoiLS2NoCD9U1bep79V9TgPOx0yZIjX9lv/GsTkyZOZN2+e21Kk7kqOuis3umzZMoKCgjwuOdoULTnqf+Xl5WzcuJHVq1fzxRdfICKMHj2aBx98kMmTJxMVFeXvEFUn57OSo/6gJUe15GhHZ7FY+Oyzz8jOzmbjxo1UVFQwYMAAsrKyyMrKqjNwQClv8EvJUdU6WnK0a8nPzyc7O5u1a9dy4sQJoqKimD59OjNmzOCqq67SAQLKLzRBdCJacrRjOX36tKMa2969ewkMDGTcuHE89dRTTJgwwWs3aSrVUpoglPKRtLQ0SktLGyyPiopi1KhRbNmyhZqaGoYPH85TTz3FtGnTiImJ8UOkSrmmCUIpH3GVHADOnz/Pnj17uO+++8jKymqXd+srBZoglPKLjz76iMDAQH+HoVSjNEEo5UVVVVV8+OGHrFy5stH1NDmojkAThFJekJuby8qVK3n33Xc5e/asVmBTnYLPEoQxpj+wFOgLWIFFIvKiMeaPwK1AFVAA3C8iZ11sXwicByxAjbtxukr5y5kzZ1i3bh2rVq1i3759BAcHk5mZye23305aWprOfKs6PF/2IGqAJ0RkpzEmCthhjNkIbASeFJEaY8x/Ak8C/+JmHzeKiOsrfUr5gcViYdu2baxatYoPPviA6upqhg8fzu9//3tuueUWevbs6Vg3NjbW5YXq2NjYtgxZqRbzWYIQkRKgxP78vDEmF4gXkfedVvsM+L6vYlDKW7755htWrVrF6tWrOX78OD169ODOO+/kjjvucHsX+bZt29o4SqW8q01KjhpjEoGrgPoFBR4ActxsJsD7xpgdxpg5jex7jjFmuzFm+8mTJ70Sb15eHk8++SR5eXle2V9TpUWfe+45Lr/8cq644gpGjhzpcd0FX5cd7eolRy9cuMCqVav40Y9+xKRJk1i0aBFDhw7lxRdfZMuWLfzud7/TKUZUp+bzi9TGmEhgJfBzESlzWv5bbKeh3nSzaZqIFBtjLgM2GmP2i8jH9VcSkUXAIrDNxdTaePPy8pg/fz5VVVXMnz+fuXPn+nSc+qeffsq6devYuXMnoaGhlJaWejyDrC/LjnbVkqMiws6dO1m5ciU5OTlcvHiRxMREfvnLXzJjxgy9+Ky6FJ8mCGNMMLbk8KaIrHJafh8wDZgobmYLFJFi++MJY8xq4FqgQYLwJufkALRJkigpKSE2NtYxrYKn56fdlR3VkqMtc/z4cdasWcPKlSspLCwkIiKCKVOmcPvtt2uZTtVl+XIUkwGWALki8men5ZOxXZQeLyIX3WzbDQiwX7voBkwCnm1tTG+//TZFRUUu37t48SJHjx5tUEmtqqqKv/zlL8THxxMREdFgu4SEBGbNmtXimCZNmsSzzz7LkCFDyMzMdBQAaoq7sqNactRztfcsrFq1ii1btmC1Whk1ahQPP/wwN998M926dfN3iEr5lS97EGnAPcBuY8wu+7KngL8CodhOGwF8JiKPGGPigMUiMhXoA6y2vx8ELBORDT6MlWPHjrksswm20w7Hjh1zfJtuLnffPo0xREZGsmPHDrZs2cLmzZuZNWsW8+bNY/bs2Y3u01XZ0WHDhmnJUQ/s37+fd955p849C3PmzOG2226rU2dDqa7Ol6OYtgKujijr3axfDEy1Pz8ENKyV2UqNfdOvf3rJWUhISKtOMzVVWjQwMJCMjAwyMjL43ve+x+uvv95ognBXdnTmzJlactSNs2fPsm7dOlauXOnyngW9s1mphvROarvU1FTmzp3bIEm0NjlA46VF8/LyCAgIcNRx2LVrFwMHDgRclxwF92VHN2zYoCVHnbi7Z+F3v/sd06ZNq3PPglKqIU0QTuonCW8kh1ruSovu2LGDxx57jLNnzxIUFERycjKLFi1yW3IU3Jcdveaaa7TkKN/ds5Cdnc2xY8c8umdBKdWQlhx1IS8vj9dee43Zs2f7bSpmLTnaPBcuXOC9995j1apVfPnllwQEBJCens4dd9zBhAkTCAkJ8Wt8SrVXjZUc1QTRicTHx3PkyBECAtrk/kcHf7Wxu3sWbr/9dr1nQSkPaU3qLqKrlBzVexaUahuaIFSH4O6ehTlz5jB58mS9Z0EpH9AEodq1/fv3s3LlStauXav3LCjVxrpEghARPe3gI764hlV7z8KqVavYu3cvwcHBTJw4kTvuuEPvWVCqDXX6BBEWFsapU6eIiYnRJOFlIsKpU6ccN+a1hsVi4ZNPPmHlypV6z4JS7USnTxAJCQkUFRXhranAVV1hYWEkJCS0eHu9Z0Gp9qvTJ4jg4GDHlBaqfXB3z8KTTz6p9ywo1Y50+gSh2gets6BUx6MJQvmUq3sWJk+ezB133KH3LCjVzmmCUF5XVVXF5s2bWblypd6zoFQH1mSCMMaMAsYCcUAFsAf4QERO+zg21cHUv2fhsssu03sWlOrA3CYIY8xs4GfAYWAHkAeEAenAvxhj9gC/F5Fv2yBO1Q6kpaVRWlraYHm3bt1ITEysc8/C7bffTnp6ut6zoFQH1lgPohuQJiIVrt40xowEUgBNEF2Eq+QAtlFJVqtV71lQqpNxmyBE5G+NbSgiu7wejeqwsrOz/R2CUsrLPJ4X2hhzqzHmc2PMLmPMT3wZlGpfKisreeONN/wdhlKqjblNEMaY+sWM7wGuA64GHm1qx8aY/saYzcaYXGPMXmPM4/blvYwxG40x+fZHl+cjjDGTjTF5xpiDxpjfeP5PUt5y4cIFlixZwsSJE/nDH/7g73CUUm2ssR7ET4wxi4wxfe2vjwDPAc8CxR7suwZ4QkSGYUssPzXGDAd+A2wSkRRgk/11HcaYQOBvwBRgOHCXfVvVBsrLy3nppZeYMGECzz//PCkpKdqDUKoLauwaxMP2XsTLxpjtwO+BG4AI4N+a2rGIlAAl9ufnjTG5QDyQBWTYV3sd+D/gX+ptfi1wUEQOARhj3rJvt8/Tf5hqvrNnz/LGG2+wdOlSysrKGD9+PI8++ihXXXUVALGxsS4vVMfGxrZ1qEqpNtDofRAi8hWQZYy5FVgLvC4izf4qaYxJBK4CPgf62JMHIlJijLnMxSbx2HostYqAMc39XOWZ06dP8/e//50333yTCxcukJmZyaOPPsqIESPqrLdt2zY/RaiU8ofG7oN4BHgYEOB5YDK2007vAX8QkS2efIAxJhJYCfxcRMo8nFrB1UouCw8YY+YAcwAGDBjgyb6V3YkTJ3j11Vd56623qKysZPLkyTzyyCMMHTrU36EppdqBxnoQPxGRK4wxIcCnIvIW8FdjzBvYTjc1mSCMMcHYksObIrLKvvi4MaafvffQDzjhYtMioL/T6wTcXPcQkUXAIoBRo0Z5v3pNJ1RSUsIrr7zCP/7xDywWC9OmTePhhx8mKSnJ36EppdqRxhLEUWPMvwHhwP7ahSJyBvhlUzs2tq7CEiBXRP7s9NZa4D5gnv1xjYvNvwRSjDGDgKPAncAPm/pM1bgjR46waNEiVq9ejYgwY8YMHn74Ye15KaVcaixBZAE3A9XAxhbsOw3b0Njdxphd9mVPYUsMK4wxD2K7C/sHAMaYOGCxiEwVkRpjzFzgPSAQeFVE9rYgBgUcOnSIl19+mXfffZeAgAB+8IMf8NBDDxEXF+fv0JRS7ZhxV1PYGJMoIoVuN7T1EOJFpMhHsTXbqFGjZPv27f4Oo904cOAAL730EuvXryc0NJRZs2bx4IMPau0FpZSDMWaHiIxy9V5jPYg/GmMCsJ0C2gGcxDZZXzJwIzAReBrb9QLVjuzbt48FCxawceNGIiIi+PGPf8z9999PTEyMv0NTSnUgjd0H8QP7zWk/Ah4A+gEXgVxgPfCciFS2SZTKI1999RULFizg//7v/4iKiuInP/kJ9957r06ep5Rqkabug9gH/LaNYlEt9OWXX7JgwQI++eQTevTowc9//nPuvvtuoqKi/B2aUqoD04pyHZSI8Omnn7JgwQK+/PJLYmJi+NWvfsVdd92lFduUUl6hCaKDERE++ugjFi5cyK5du7jssst46qmnmDlzJuHh4f4OTynViWiC6CCsViubNm1i4cKF7N27l/j4eJ555hnuuOMOQkJC/B2eUqoT8qQm9UrgVSBHRKy+D0k5s1gsbNiwgZdeeokDBw4wcOBA/v3f/53p06cTHBzs7/CUUp2YJz2IhcD92KbZ+Afwmojsb2Ib1Uo1NTWsW7eOl156icOHD5OUlMQf//hHpk6dSlCQdvyUUr7X5JFGRD4APjDGRAN3ARuNMUeAV4D/EZFqH8fYpVRVVbFmzRpefvlljhw5QmpqKv/1X//FzTffTECAxwUAlVKq1Tz6KmqMiQHuxjZ1xj+BN4F0bHMpZfgquK7k0qVLvPPOO7zyyiuUlJQwYsQInnzySW688UZNDEopv/DkGsQqYCjwBnBrbS0H4G17ISHVChcvXuTtt99myZIlnDx5kquuuopnn32WsWPH4uHU6Eop5RON1YPoZ08G80XkQ1fruJu/QzWtvLycZcuW8fe//53Tp08zZswYXnjhBcaMGaOJoZPJy8vjtddeY/bs2aSmpvo7nHZP26v9aGyyvhygJ7aSoBuArSJS03ahNV9HmKyvrKyMpUuXsnTpUs6dO8fYsWN59NFHueaaa/wdmvKBvLw85s+fT1VVFSEhIcydO1cPeo3Q9mp7jU3W5zZB2DcMw3aNYQq26bu/xZYsNojIt94PtXXac4I4ffo0r7/+Ov/zP/9DeXk5EyZM4NFHH+WKK67wd2jKR5wPdrX0oOeetpd/tDhBuNjRIGzJYjLQV0Su9U6I3tEeE8TJkycdZT0rKiqYNGkSP/nJT7SsZydTXV1NRUUFlZWVVFZWkpeXR3Z2NjU1DTvdgYGBXH/99fTu3RsRcfyA7YZIwOVr5/Uae+3JOi3dxlfxVVZWUlpa6njPmTGG5ORkYmJiCA4OJjg4mJCQkDrPa1/XX17//eDgYIKCgjrVadzWnpLzSoIwxnSn7jWLchGpcre+P7SnBHHs2DEWL17MihUrqK6u5pZbbuGRRx4hOTnZ36EpOxGpc2B3PsDXX9bYOpWVlS4TQWvVHsSMMY6fpl57uk3tyLiW7tvTfXgaz/79+6mudj9iPiAggJ49e1JdXU1VVRXV1dVYLJYWt6ur5FI/qbh6bCzxuNvGlwnJG6fkWloPonbjh4FngQqgNpuIiAxuVhRdxJEjR3jllVdYtWoVIkJWVhZz5swhMTHR36F5jb8vIooIly5dcnswr33u7mDuvH7tN+DGBAUFER4eTlhYGGFhYYSHh9OzZ0/69etHeHh4g/fCwsI4efKk2x5EcHAwjzzyCEOGDGn0YNqVuDq9VMvdgc9isdRJGM7P6y9rarnzexUVFS735a2E5EkCaizx1D4vKSnhH//4hyOxVlVVMX/+fK+ekmuyB2GMyQeuF5FSr3yiD/mzB1FYWMjLL7/MmjVrCAgI4I477uChhx4iISHBL/H4Smu+sVitVseBvTXf1isqKlyeiqgvODjY5cHb3WP9ZeHh4YSGhrZ4ShM9p9487b29rFZrg6TR3ITUnG1bmpCa22atOsVkjNkA3C4iF5sdaRvzR4I4ePAgCxcuZP369QQHBzNz5kx+/OMf07dv3zaNoy3UVqpzPhUQGBhIWloa3bt3b/IAf+nSJY8O7KGhoQ0O2M4HbXcH/PrrBAYG+rI5PKKjcppH2+s7zgmpfvJYuHAh58+fd7ttr169+I//+A+PPqe1CeIq4O/A58Cl2uUi8jOPPr0NtWWCyM3NZeHChbz//vuEh4dz1113cf/999O7d+82+Xxvqb1AeO7cOc6dO0dZWZnjufPr06dPU1nZeAHB5hzUXX1bDwsLIzQ0tF0c2L3J36fkOhptr6a15JScO61NEF8AW4HdgOOErYi83sR2rwLTgBMiMsK+7G2gNuoewFkRGeli20LgPGABajy9Ia8tEsTXX3/NggUL2Lx5M5GRkdx9993cd9999OrVy6ef21xWq5Xy8vI6B3tXSaCsrMzlL1lQUBDR0dF0796d6OhocnNzuXTpkotPsmnONxalVOt565Rcqy5SYztA/9LjT/vOa8B8YGntAhGZ5RTUn4BzjWx/Y3u67rF9+3YWLlzI1q1biY6O5mc/+xn33HMP3bt3b9M4qqurGz3Y1z4/f/68ywuw4eHhREdHEx0dzaBBgxzPnZNBdHQ0ERERdS6UNvWNZfbs2b78Zyul6klNTWXu3Lk+PSXnSQ/iOeAb4F3qnmI63eTOjUkE1tX2IJyWG2w33U0QkXwX2xUCo5qbILzdgxARPvvsMxYsWMAXX3xBr169uP/++/nhD39IZGSkVz+noqLC7Wke52UXLza8FGSMISoqyu3B3nlZa4oLtfeLiEp1RX69D8IYc9jFYo+GuTaSIMYBf3YblO0zz2AbVvuyiCxq5DPmAHMABgwYcM0333zTVFgN1G9gEeHjjz9m4cKF/POf/6R37978+Mc/ZubMmURERHi8X4vF4vFpHldjwIODg90e9J2fR0VFtdmMr3oRUanOxWt3UrfggxNxnSAWAgdF5E9utosTkWJjzGXARuAxEfm4qc9rbg8iLS2Nmpoavve97xEYGIjFYmH37t2cP38ei8VCv379eOihh/j+979PaGioY7uqqiqPT/O4at+IiIgmv+1HR0cTFhbWLsfD60VEpTqP1vYgwoCfYKv/IMAW4CURaXxIC64ThDEmCDgKXCMiRR7s4xlsd22/0NS6zU0QY8aMcSSHWhaLhcOHD5OVlUVycjLl5eUNTvm4Gs0TEBBA9+7d3R7sa5d3795dS4UqpdqN1l6kXoptRNF/21/fha02xA9aGE8msN9dcjDGdAMCROS8/fkkbHdye1VeXl6D5AC2cf3Jycns3buXvXv3EhIS4jjIJyQkMHz4cJcH/sjISC3so5TqVDxJEKkicqXT683GmK+a2sgYsxzbTLCxxpgi4GkRWQLcCSyvt24csFhEpgJ9gNX2UytBwDIR2eDJP6Y5XnvttUbH20dHR/Pss88SFhbm7Y9WSqkOwZME8U9jzHUi8hmAMWYMsK2pjUTkLjfLZ7tYVgxMtT8/BFxZfx1vmz17Nn/84x9dJgmLxcKDDz6oyUEp1aV5ck5kDPCJMabQPvz0U2C8MWa3MeZrn0bnQ6mpqezevbvBfCe1F6r14qtSqqvzpAcx2edR+ElQUBC7d+9uMIopKMiTZlFKqc6tySOhiDT/xoIOYts225kyHbaplFIN+fQ+iLbWngoGKaVUR9DaYa5KKaXaoa1bt7JixQpKS0uJjY1l5syZpKene23/miCUUqoD2rp1K4sXL3bMjVZaWsrixYsBvJYk9M4upZTqQM6ePcsXX3zBq6++2mB25aqqKlasWOG1z9IehFJKtVMiQnFxMQcOHCAvL4+8vDyOHz/e6DanTp3y2udrglBKqXaipqaGw4cPO5LBgQMHHKVFIyMjSU1NZeLEiaSmpvLXv/7VZTKIiYnxWjyaIJRSyk8uXLhAfn6+IyEUFBQ4pv7v06cPV111FampqaSmptKvX786szvPmjWrzjUIsNVnmTlzptfi0wShlFJtpLS01JEM8vLyKCoqQkQICAggMTGRzMxMR0KIjo5udF+1F6JXrFjBqVOniImJ0VFMSinVEVitVo4cOVLndFHtKaGwsDBSUlIYM2YMqampJCUltWjut/T0dK8mhPo0QSillBdcunSJgoICxwXlAwcOUFFRAUDPnj1JTU1l2rRppKam0r9//0Znk24vNEEopVQLlJWVORJBXl4ehw8fdkz+2b9/f2644QbH6aLY2Nh2WR2yKZoglFKqCSLCsWPH6iSEkpISwDbpZ1JSErfccgupqamkpKQQGRnp54i9QxOEUkrVU1NTQ2FhYZ2EUFZWBtiGmw4ZMoSMjAxSU1MZNGhQpy0jrAlCKdXlXbx40THc9MCBAxw8eNAxfPSyyy7jyiuvJDU1lSFDhhAXF9dlygtrglBKdTmnTp2qM7ro22+/RUQwxpCYmMiECRMYMmQIqamp9OzZ09/h+o0mCKVUp2a1WikqKqpzuqi0tBSA0NBQkpOTue2220hNTSU5OZnw8HA/R9x++CxBGGNeBaYBJ0RkhH3ZM8BDwEn7ak+JyHoX204GXgQCgcUiMs9XcSqlOpeqqqoGw00vXrwIQI8ePRgyZAhTpkwhNTWVgQMHdojhpv7iyx7Ea8B8YGm95X8RkRfcbWSMCQT+BtwEFAFfGmPWisg+XwWqlOq4ysrK6kxXcejQIcdw0/j4eMfNaKmpqVx22WUdcripv/gsQYjIx8aYxBZsei1wUEQOARhj3gKyAE0QSnVyTRXAERFOnDhRZ7qK4uJiAAIDA0lKSmLq1KkMGTKEIUOGEBUV5a9/Sqfgj2sQc40x9wLbgSdE5Ey99+OBI06vi4Ax7nZmjJkDzAEYMGCAl0NVSrUVdwVwTpw4QXh4uCMhnDt3DoCIiAhSU1MZN24cQ4YMYfDgwYSEhPjzn9DptHWCWAj8GyD2xz8BD9Rbx1X/z23hbBFZBCwCW01q74SplGprK1ascFkA55133gGgd+/ejBgxwnG6KD4+vssMN/WXNk0QIuKodGGMeQVY52K1IqC/0+sEoNjHoSml2pjVaqW4uJj8/Hzy8/MdI4tc+e///m+v1jlQnmnTBGGM6SciJfaXtwF7XKz2JZBijBkEHAXuBH7YRiEqpXzkwoULHDx40JEQCgoKHKOLunXrRnBwsKMWgrPY2FhNDn7iy2Guy4EMINYYUwQ8DWQYY0ZiO2VUCDxsXzcO23DWqSJSY4yZC7yHbZjrqyKy11dxKqW8z2q1cvTo0ToJ4ejRowAYY0hISOC6664jJSWFlJQU+vbtyyeffOLzAjiqeYxI5zltP2rUKNm+fbu/w1CqyykvL+fgwYOOhHDw4EHHVNeRkZEkJyc7kkFSUpLbm9FqRzH5qgCOasgYs0NERrl8TxOEUqo5au9Mdu4d1A41NcYwYMCAOgmhb9++eu9BO9ZYgtCpNpRSjSovL3f0CmqvHTj3DlJSUkhPTyclJYXBgwfrVBWdiCYIpZRDbe+gtmeQn5/vqHtQ2ztIS0tz9BC0d9C5aYJQqgs7f/58nVNFhw4dcvQOoqKiSElJYdy4cSQnJ7e4brLquDRBKNVFWK1Wjhw5Uud0UW3vICAggAEDBjhOFSUnJ9OnTx/tHXRxmiCU6qTKysrqjCoqKCigsrISgO7du5OSksL48eNJSUlh0KBB2jtQDWiCUKoTsFgsjt5BbUI4duwY8F3vYOzYsY6RRTqrqfKEJgilOqDaKa6deweXLl0CvusdZGRkOEYWhYaG+jli1RFpglCqnavfO8jPz+f4cdu0ZoGBgQwYMMBxqiglJYXevXtr70B5hSYIpdqZc+fONbh2UNs7iI6OJiUlhQkTJjiuHWjvQPmKJgilfKipAjgWi4Vvv/22Tu/gxIkTgK13MHDgQMepopSUFGJjY7V3oNqMJgilfMRdAZzDhw8TFBREfn4+hw8fdvQOevToQUpKChMnTnRcO9ACOMqfNEEo5SNvv/22ywI4OTk5BAYGkpiYyI033ui4K1l7B6q90QShlBdYLBaOHj1KQUEBhw4doqCggFOnTrldf8mSJdo7UO2eJgilmklEOH78uCMRHDp0iMLCQsepooiICMekdbXTVjiLjY3V5KA6BE0QSjXhzJkzdXoGhw4d4sKFCwAEBwczaNAgbrzxRgYPHkxSUhJ9+vQhICCgwTUI0AI4qmPRBKGUk/Lycg4fPuxIBIcOHeL06dOA7Y7k/v37c+2115KUlMTgwYNJSEggKMj1n1HtaCUtgKM6Kk0Qqsu6dOkShYWFdXoGtdNTAPTt25dhw4Y5egYDBw5s9j0H6enpmhBUh6UJQnUJNTU1FBUV1ekZHDlyBKvVCkCvXr0YPHgw48ePZ/DgwQwaNIjIyEg/R62Uf/ksQRhjXgWmASdEZIR92R+BW4EqoAC4X0TOuti2EDgPWIAad+XwlHLFarVy7NixBheRq6urAejWrRuDBw9m+vTpjlNFPXv29HPUSrU/vuxBvAbMB5Y6LdsIPCkiNcaY/wSeBP7FzfY3ikipD+NTnYCIcPr06ToXkQ8fPszFixcBCA0NJTExkczMTJKSkkhKStKZTJXykM8ShIh8bIxJrLfsfaeXnwHf99Xnq87p/PnzdXoGhw4d4uzZs4Btaor+/ftz/fXXO64bxMfHExgY6N+gleqg/HkN4gHgbTfvCfC+MUaAl0VkkbudGGPmAHMABgwY4PUglf9UVlZSWFhIQUGBIyHUzlMEEBcXx4gRIxyniQYOHKj3FyjlRX5JEMaY3wI1wJtuVkkTkWJjzGXARmPMfhH52NWK9uSxCGDUqFHik4CVz9XU1PDtt9/W6R0UFRUhYvsvjY2NZfDgwUyYMIGkpCQGDRpERESEn6NWqnNr8wRhjLkP28XriVL711+PiBTbH08YY1YD1wIuE4TqeKxWKyUlJXWuG3z77beOi8iRkZEkJSUxevRox6mi6OhoP0etVNfTpgnCGDMZ20Xp8SJy0c063YAAETlvfz4JeLYNw1ReJCKUlpY6rhfUXkSunYIiNDSUwYMHM2nSJMepIi14o1T74MthrsuBDCDWGFMEPI1t1FIottNGAJ+JyCPGmDhgsYhMBfoAq+3vBwHLRGSDr+JUzdNUfYOysrIG01KUlZUB39U3SE9Pd/QM4uLiCAgI8Nc/RynVCOPmLE+HNGrUKNm+fbu/w+i0XM0tFBwczKhRo7BarRQUFFBaahuZbIwhLi7OMbR08ODBDBgwgODgYH+Fr5RywRizw929ZnontfLY8uXLG9Q3qK6u5tNPP6V3794kJSVx0003kZycTGJiIuHh4X6KVCnlDZoglEsiwsmTJ9m3bx/79+9n3759nDlzxu36L774YhtGp5RqC5ogFPBdjYPc3FzHT23Bm6ioKIYOHUpFRYVjmmtnsbGxbR2uUqoNaILookSEY8eOsW/fPnJzc9m/f79jWuvu3bszdOhQbr31VoYNG0Z8fLzWN1CqC9IE0UWICMXFxXV6CLVTVERHRzNs2DCGDRvG8OHDiYuLcznMVOsbKNW1aILopESEo0ePkpub67iOcO7cOQB69uzJ8OHDHUmhX79+Ht93oPUNlOo6NEF0ElarlaKiIkfvYP/+/Y77D3r16sWIESMYPnw4Q4cOpW/fvnojmlKqSZogOiir1cqRI0ccvYPc3FzKy8sB20XjK6+80tFD0OmtlVItoQmig7BarXzzzTd1egi1I4p69+7N1Vdf7Tht1Lt3bz9Hq5TqDDRBtFMWi8WREPbt20deXp6jCE6fPn0YPXq0o4egw0yVUr6gCaKdsFgsHD582NFDyMvLc0xo17dvX8aMGeNICDExMX6OVinVFWiC8JOamhoOHTrkuH6Ql5dHZWUlYCuEc8MNNzgSgtZL7tguXbrEvn37GD58OKGhof4Op93T9mo/NEG0kZqaGgoKChw9hAMHDnDp0iUA4uPjSU9Pd4wy6tGjh3+DbYL+AduGETv/WK3WBstqfwoLCzl37hz5+fkMHDgQoM6ggeY8b+l2nuyjvahtr8LCQlJTU/0dTpemCcJHqqurKSgocNypnJ+f77gDuX///owfP55hw4YxdOjQDlcMp6k/YHcHSk8Opk293162bYnS0lLHbLftUVskqsberx2FV6ukpISSkhLANjLPGENAQAABAQGO5+4eW/qe82N7TJ6u+PILmyYIL6mqqmqQEGorpA0cOJAbb7zRkRC6d+/u52i/IyJYLBZqamocP9XV1XVe1y5zrgcNdf+AAwMDW30Aba3aP+raH+c/9Mbec7Wep9s29n5NTQ0nTpygrKwMEcEYQ/fu3enTp49j2nPntmrqeXPW9cZ2bf3ZwcHBXLhwoc5ULkFBQYSEhFBRUeFI0s6Ptc99oTUJprUJqjlJzJc9Lk0QLVRVVUV+fr5jlFFBQQHV1dUYYxg4cCCZmZmOhBAZGenTWJwP8vUP7u4O9s6vm/oDCwoKIigoiIiICKqrqx2JDyAsLIzo6GiCg4NbdED19GDryXvtUXl5OefOnSMgIACr1UpERARxcXH+DqvdysvLo6SkxNFevXv3bvKg56qn55w8Gnts6XvO61gsFo+297XaL2wBAQGMGzfOK/vUBOGhyspKR0LIzc2loKCAmpoajDEkJiYyadIkhg4dSmpqaosSgog0eSBv7HVTgoKCCA4Odhzsw8LCGixz9TooKKjOwbf+H3DPnj31PHEjqquriYuLIy4ujuLi4gb1NFRdLWkv5y8IgYGBvg6xRXyRxKqqqjhz5oxj+HtAQACxsbEkJSV5Le4unyC2bt3K2rVrGTlyJP/85z/JysoiPT2dyspKDhw44LhTuaCgAIvFQkBAAImJiUyePJlhw4aRmppKREQE8N1B/uLFi80+2Dd1kDfGNDiQh4eHN3pgr10WGBjotW/YesBrnhEjRjieDxkyxI+RdAydtb18lcRq74+q/cIWGBjo1esQXbrkaO301aNHj2b48OHs27ePzz//nL59+3LhwgWCgoIIDw+nf//+JCQk0KdPH8cII1cHe4vF0ujn1R7km/rW7upg782DvFKqc9izZw8hISF1vrA5J1lPmEZKjnbpBLFp06ZmZ/OAgACPvrW7et1eu79Kqa6rsQThs1NMxphXgWnACREZYV/WC3gbSAQKgZki0qCOpTFmMvAiEAgsFpF5vohx2bJlXHfddSQlJTm6aKdPn+bQoUPce++9Lg/2AQEBvghFKaXaHV8e7V4DJtdb9htgk4ikAJvsr+swxgQCfwOmAMOBu4wxw30RYLdu3Rwjj2ovOB8/fpyioiL69u1LTEwM0dHRdOvWjdDQUE0OSqkuxWdHPBH5GDhdb3EW8Lr9+evADBebXgscFJFDIlIFvGXfzutmzpxJt27d2LdvH9nZ2ezbt4/IyEgtoamUUvi2B+FKHxEpAbA/XuZinXjgiNPrIvsyl4wxc4wx240x20+ePNmsYNLT00lKSiIvL48zZ86Ql5dHUlKSVkxTSina5zBXV0N13F5JF5FFwCKwXaRu7odpCU2llHKtrXsQx40x/QDsjydcrFME9Hd6nQAUt0FsSimlnLR1glgL3Gd/fh+wxsU6XwIpxphBxpgQ4E77dkoppdqQzxKEMWY58CmQaowpMsY8CMwDbjLG5AM32V9jjIkzxqwHEJEaYC7wHpALrBCRvb6KUymllGs+uwYhIne5eWuii3WLgalOr9cD630UmlJKKQ/owH6llFIudaqpNowxJ4FvWrh5LNB+q7m0P9pezaPt1TzaXs3TmvYaKCK9Xb3RqRJEaxhjtrubj0Q1pO3VPNpezaPt1Ty+ai89xaSUUsolTRBKKaVc0gTxnUX+DqCD0fZqHm2v5tH2ah6ftJdeg1BKKeWS9iCUUkq5pAlCKaWUS502QRhj+htjNhtjco0xe40xj9uX9zLGbDTG5Nsfe9qXx9jXLzfGzHfaT5QxZpfTT6kx5r/89M/ymRa017VObfKVMeY2p31dY4zZbYw5aIz5q+mExbSb215O2w2w/479P6dl2l4Nf78SjTEVTr9jLzntS9vLxe+XMeYKY8yn9vV3G2PC7Mtb3l4i0il/gH7A1fbnUcABbBXqngd+Y1/+G+A/7c+7AenAI8D8Rva7Axjn739fO2ivCCDIadsTTq+/AK7HNnV7DjDF3/8+f7eX03YrgX8A/89pmbZXw9+vRGCPm31pezVsryDga+BK++sYILC17dVpexAiUiIiO+3Pz2Ob+C8eN1XtROSCiGwFKt3t0xiTgq3I0RbfRe4fLWivi2KbWBEgDHvNDvs07t1F5FOx/XYuxXXlwA6tue0FYIyZARwC9jot0/Zqusqkg7aX2/aaBHwtIl/ZtzklIpbWtlenTRDOjDGJwFXA53hW1c6du4C37Q3daXnaXsaYMcaYvcBu4BF7wojHVtOjVqMVATsDT9rLGNMN+BfgX+ttru3l/u9xkDHmn8aYj4wxY+3LtL1ct9cQQIwx7xljdhpjfm1f3qr2ao8V5bzKGBOJrVv/cxEpa+XpyjuBe7wSWDvVnPYSkc+By40xw4DXjTE5NLMiYEfXjPb6V+AvIlJebx1tL9dKgAEicsoYcw2QbYy5HG0vd6sGYTtFPhq4CGwyxuwAylys63F7deoEYYwJxta4b4rIKvvi48aYfiJSYtxXtXO1ryuxnWPf4aNw/a6l7SUiucaYC8AIbN9QEpze7rQVAZvZXmOA7xtjngd6AFZjTKV9e22veu0lIpeAS/bnO4wxBdi+Jevvl+vfryLgIxEptW+7Hrga+B9a0V6d9hST/Ur9EiBXRP7s9JYnVe1cuQtY7r0I25fmtpexVfwLsj8fCKQChfZu73ljzHX2fd6L523cYTS3vURkrIgkikgi8F/Av4vIfG0vt79fvY0xgfbng4EU4JC2l9vj13vAFcaYCPvf5XhgX6vbq62vzrfVD7bulmC7sr/L/jMV29X9TUC+/bGX0zaFwGmgHFtGHu703iFgqL//Xe2lvbCdattrX28nMMNpX6OAPUABMB/7Hfud6aclv19O2z5D3VFM2l4Nf7/usP9+fWX//bpV26vJ49fd9jbbAzzvjfbSqTaUUkq51GlPMSmllGodTRBKKaVc0gShlFLKJU0QSimlXNIEoZRSyiVNEEoppVzSBKGUUsql/w+mmKDgXgGTtwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pmoysf       = np.zeros((2,4))\n",
    "pmoys0f      = np.zeros((2,4))\n",
    "pmoysf[0,:]  = yy_eu\n",
    "pmoysf[1,:]  = yy_us\n",
    "pmoys0f[0,:] = yy_eu0\n",
    "pmoys0f[1,:] = yy_us0\n",
    "\n",
    "\n",
    "xx      = np.zeros((2,5))\n",
    "xx[0,0] = PMoY10[0]\n",
    "xx[1,0] = PMoY10[-1]\n",
    "zz      = np.zeros((2,5))\n",
    "zz[0,0] = PMoY10[0]\n",
    "zz[1,0] = PMoY10[-1]\n",
    "for jj in range(2):\n",
    "    for ii in range(horizon):\n",
    "        xx[jj,ii+1] = pmoysf[jj,ii] \n",
    "        zz[jj,ii+1] = pmoys0f[jj,ii]\n",
    "    \n",
    "lyhorizon = ('2017','2030','2040','2050','2060')\n",
    "\n",
    "plt.plot(np.arange(len(lyhorizon)),list(xx[0,:]),'o-',label=r'EU, $\\Delta p \\neq 0$',color = '0.35')\n",
    "plt.plot(np.arange(len(lyhorizon)),list(xx[1,:]),'s-',label=r'US, $\\Delta p \\neq 0$',color = '0.15')\n",
    "plt.plot(np.arange(len(lyhorizon)),list(zz[0,:]),'*-',label=r'EU, $\\Delta p = 0$',color = '0.75')\n",
    "plt.plot(np.arange(len(lyhorizon)),list(zz[1,:]),'D-',label=r'US, $\\Delta p = 0$',color = '0.4')\n",
    "plt.xticks(np.arange(len(lyhorizon)), labels=lyhorizon)\n",
    "plt.legend() #(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.ylabel('pm/y (%)')\n",
    "plt.savefig('forecasts.eps',dpi=600)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b73ef8d1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
